#include "cond.h"
#include <iostream>

#define NE    2	    /* ne equations for ne adjustable dependent variables */
#define NB    1	    /* nb boundary conditions */

inline double RLX(double x, double delta) { return (delta>0)? exp(- (x/(1e-3*delta)) ) : 0;}


/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

class Device {

 protected:

  void q_eraser();
  void set_strain(short flag = 1);

  /* Schrodinger related */ 

  int find_hetero_cb( int **, int **, FILE*);
  int find_hetero_vb( int **, int **,FILE *);
  Step * step_approximation_cb( int &,int,int,int,int,FILE *,FILE *);
  Step * step_approximation_vb( int &,int,int,int,int,FILE *,FILE * );
  double ** transmission_coefficient( int,Step *,double,double,double,int &,FILE *,int left = 0,int right = 0	);
  double ** transmission_coefficient_bidir( int,Step *,double,double,double,int &,FILE *,int left = 0,int right = 0	);
  Eigenvalue ** eigenvalues(int,Step *,FILE *,FILE *);
  int hole_wavefunction_normalization(double ,FILE *);
  int electron_wavefunction_normalization(double ,FILE *);
  void wavefunction_n( Eigenvalue ** eigen_store,int step_num,Step * band,FILE * fp, char str1[30] );
  void wavefunction_p( Eigenvalue ** eigen_store,int step_num,Step * band,FILE * fp, char str1[30] );
  void Schrodinger(short);


  /* grid interface */

  int LayerToMesh(int k, int i);
  void MeshToLayer(int j, int & k, int & i);
  void incr(int & k,int & i) {
    if(k>=num_lyr) nrerror("DEVICE:INCR. Out of Range!");
    if(i<lyr[k].numx()-1) i++;
    else { i=0; k++; }
  }
  void decr(int & k,int & i) {
    if(k<0 || (k==0 && i==0)) nrerror("DEVICE:DECR. Out of Range!");
    if(i>0) i--;
    else { k--; i=lyr[k].numx()-1;}
  }
  int test_upper_limit(int k,int i,int k_limit=-1, int i_limit=-1) {
    if(k_limit==-1 && i_limit==-1 ) {
      k_limit = num_lyr -1;
      i_limit = lyr[num_lyr-1].numx()-1;
    }
    if(k>k_limit || (k==k_limit && i>i_limit)) return 0;
    else return 1;
  }
  int test_lower_limit(int k,int i,int k_limit=-1, int i_limit=-1) {
    if(k_limit==-1 && i_limit==-1 ) {
      k_limit = 0;
      i_limit = 0;
    }
    if(k<k_limit || (k==k_limit && i<i_limit)) return 0;
    else return 1;
  }

  double phi(int);
  double phi2D(int,int);
  double charge(int,int,double);
  double dcharge(int,int,double);
  double Polar(int,int);
  double eps(int,int);
  void set_phi(int,double);
  void set_phi(int,int,double);
  double phin(int);
  double phin(int,int);
  void set_phin(int,double);
  void set_phin(int,int,double);
  double phip(int);
  double phip(int,int);
  void set_phip(int,double);
  void set_phip(int,int,double);
  double D(int);
  void set_D(int,double);
  double Jn(int);
  void set_Jn(int,double);
  double Jp(int);
  void set_Jp(int,double);
  double doping(int);
  double chi(int,int);
  double n(int);
  double p(int);
  double ni(int);
  double n1(int);
  double p1(int);
  double pn1(int);
  double np1(int);
  double chi(int);
  double Eg(int);
  double Nc(int);
  double Nv(int);
  double mobn(int);
  double mobp(int);
  double mobn(int,int);
  double mobp(int,int);
  double me(int);
  double mh(int);
  double nf(int,int,double);
  double dnf(int,int,double);
  double pf(int,int,double);
  double dpf(int,int,double);
  double Dn(int);
  double Dp(int);
  double Tn(int);
  double Tp(int);
  double Tntrap(int);
  double Tptrap(int);
  double SHR_Rec(int);
  double Trap_Rec_e(int);
  double Trap_Rec_h(int);
  double Opt_Rec(int);
  double QW_Rec(int);
  double dx(int);
  double dx(int,int);
  double dy(int);
  double dy(int,int);
  double dt(int);
  double dt(int,int);

  /* Transport related */

  void Poisson_solvde(short);
  void difeq( int k,int k1,int k2,int jsf,int is1,int isf,int indexv[],int ne,double **s, double **y );
  double ** ElectronDiscretization1D(short type, int & n);
  double ** HoleDiscretization1D(short type, int & n);
  double ** ElectronDiscretization1D_2(short type, int & n);
  double ** ElectronDiscretization1D_3(short type, int & n);
  double ** HoleDiscretization1D_2(short type, int & n);
  double ElectronSolver(short,double,short);
  double HoleSolver(short,double,short);
  void Current(short,short);
  double GrabCurrent(short debug);
  void Poisson_2D();
  void Continuity_2D();

 public:


  int num_lyr;
  int xnum_mesh_points;
  int ynum_mesh_points;
  int tnum_mesh_points;
  Layer * lyr;
  Condition * sim;
  Device(double PZ=-1,double PSP=-1);
  ~Device() {delete [] lyr; delete sim;}
  void dump_files(char);
  void mesh_modifier();
  void ShowMatrix(double **, int, int);
  void ShowVector(double *, int);
  int simulator();
  double field_in_qw();
  double * bias_distribution(double, short);
  double ** optical_transitions(int & num_transition);
  friend void solvde( int,double,double,double scalv[],int indexv[],int,int,int,double **,double ***,double **,Device &);

};

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf, int ic1, int jc1, int jcf, int kc, double ***c, double **s);
void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k, double ***c, double **s);
void bksub(int ne, int nb, int jf, int k1, int k2, double ***c);
void solvde( int itmax,double conv,double slowc,double scalv[],int indexv[],int ne,int nb,int m,double **y,double ***c,double **s,Device & dev  );

